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I. INTRODUCTION 

The experimental progress in controlling quantum op- 
tical and atomic systems, which has been achieved over 
the last few years, prompted ideas for new realizations 
of strongly correlated many body systems, such as ultra- 
cold gases of atoms trapped in optical latticesi"— or light- 
matter systems The latter consist of photons, which 
interact with atoms or atomiclike structures. Normally, 
the interaction between photons and atoms is very weak, 
since the interaction time is small. However, a strong 
interaction can be achieved when photons are confined 
within optical cavities. In this case, the coupling be- 
tween photons and atoms leads to an effective repulsion 
between photons, which means that it costs energy to 
add additional photons to the cavity. The arrangement 
of such cavities on a lattice, see Fig.[TJ allows the photons 
to "hop" between neighboring sites, provided the cavities 
are coupled. Quantum mechanically the coupling of ad- 
jacent cavities means that their photonic wave functions 
overlap. Due to the strong interaction between photons 
and atoms, and the introduction of a lattice of coupled 
cavities, a strongly correlated phase emerges where pho- 
tons are present. The light-matter models share some 
basic properties with the Bose-Hubbard (BH) model, ^ 
such as the quantum phase transition from a Mott phase, 
where particles are localized on the lattice sites, to a 
superfluid phase, where particles are delocalized on the 
whole lattice!^ Yet the physics of the light-matter mod- 
els is far richer because two distinct particles, namely. 
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FIG. 1: Cavities forming a one-dimensional chain lattice. The 
blue dots represent atomic systems, whereas the red wavy 
arrows indicate photons. 



photons and atomiclike excitations, are present. 

A major advantage of these man-made realizations of 
strongly correlated many-body systems is that they can 
be tailored to correspond to a many-body model, whose 
parameters can be directly controlled in the experiment. 
Furthermore local quantities, such as the particle den- 
sity at a specific lattice site, can be addressed individu- 
ally due to the mesoscopic scale of the cavities and both 
lattice size and geometry can be controlled in the fab- 
rication process. An experimental realization of these 
light-matter systems is still missing but there are sev- 
eral promising approaches, such as photonic crystal cavi- 
ties or toroidal and disk-shaped cavities If light-matter 
systems can be realized, they will undoubtedly provide 
fascinating insight in the physics of strongly correlated 
many-body systems. The realizations might be used as 
quantum simulators for other quantum mechanical prob- 
lems or even more intriguing for quantum information 
processing applications.^ 

Recently, there has been a lot of research activity in the 
field of light-matter systems. Most of the work has been 
devoted to investigate the quantum phase transition from 
the Mott to the superfiuid phase. Some basic character- 
istics of the quantum phase transition have been evalu- 
ated from small systems of a few cavities by means of ex- 
act diagonalization.i^— Results are available at mean- 
field level^j^ — as well or more accurately from analyt- 
ical strong coupling perturbation theory calculations^ 
and from simulations based on the density matrix renor- 
malization groupiiii^ (DMRG), the variational cluster 
approacbi^ and Quantum Monte Carlo Spectral prop- 
erties of light-matter systems have been investigated in 
Refs. [ii,[l9|, andllH. 

In the present paper, we study in detail the spectral 
properties of a one-dimensional light-matter system. In 
particular, we evaluate both photonic as well as atomic- 
excitation spectral functions. The investigation of both 
spectral functions allows us to characterize the polariton 
excitations in light-matter models. In addition to the 
spectral functions, we present densities of states, momen- 
tum distributions and spatial correlation functions. For 
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completeness we also show the first two lobes delimiting 
the Mott transition. 

This paper is organized as follows: in Sec. [Til we intro- 
duce the light-matter model. Section [TTll contains both 
the description of the numerical method as well as the 
exploration of the polaritonic properties. Section |IVl is 
devoted to spectral properties of the light-matter system. 
Here, we present our results for spectral functions, den- 
sities of states, momentum distributions and spatial cor- 
relation functions. Finally, we summarize and conclude 
our findings in Sec. |Vl 

II. MODEL 

From the great variety of possible theoretical descrip- 
tions of light-matter systems^ —li^i^i^ we concentrate on 
the simplest one, which consists of an array of cavities 
each of which contains a two- level system.^ The physics of 
the i-th cavity can be described by the Jaynes-Cummings 
(JC) Hamiltonian,^^ which for = 1 is given by 

H-^ = al a, + e (7+ a' + g (a, (7+ + a] a'^ , (1) 

where ujc is the resonance frequency of the cavity, i.e., 
the frequency of the confined photons, e is the energy 
spacing of the two-level system, and g is the atom-field 
coupling constant. The operator a] creates a photon 
with frequency cjc, whereas annihilates one. The two- 
level system can be mathematically described by Pauli 
spin algebra. Thus, we identify the ground state of the 
two- level system with ||^) and the excited state with 
Iti). With that the atomic raising operator is defined 
as = Iti) and the atomic lowering operator as 
= (til, respectively. In order to obtain the JC 
Hamiltonian the rotating wave approximation, which is 
justified for \ujc — e| <C c^c, has been assumed. The 
deviation between the resonance frequency and the en- 
ergy spacing of the two-level system, A = — e, is 
termed detuning. For the JC Hamiltonian the particle 
number hi = a| + a~ is a conserved quantity, as 
[Hf^^ hi] = 0. This is a consequence of the rotating 
wave approximation*^^ 

The full model consists of an array of TV cavities, which 
form a lattice and hence we refer to this model as the 
Jaynes-Cummings lattice (JCL) model. Due to the cou- 
pling of the cavities, photons are allowed to hop between 
neighboring lattice sites. This leads to the JCL Hamil- 
tonian 

^JCL = -tY, 4 a, +Y^HfO-i^Np, (2) 

where t is the hopping strength and fi the chemical poten- 
tial, which controls the total particle number Np of the 
system. The first sum with the angle brackets around the 
summation indices is restricted to nearest-neighbor sites. 
In the case of the JCL model, the particle number of a 



specific cavity hi is not conserved anymore. However, the 
total particle number Np = hi is a conserved quan- 
tity. In summary, the JCL Hamiltonian can be rewritten 
as 

+ g "^{^i + 4 err) - (/i - cjc) Np . (3) 

i 

From Eq. (|3|) and from the fact that we consider the cou- 
pling strength g as unit of energy, it follows that the 
physics only depends on three independent parameters, 
namely the hopping strength the detuning A and the 
modified chemical potential fi — Uc. In order to fulfill the 
condition for the rotating wave approximation the reso- 
nance frequency uJc has to be large in comparison with the 
detuning A, which can be always satisfied theoretically as 
solely the difference between the chemical potential and 
the resonance frequency appears in the grand-canonical 
Hamiltonian H'^^^. 



III. METHOD 

In order to investigate the properties of the JCL model, 
we use the variational cluster approaches (VCA), which 
has been formulated for bosonic systems in Ref. [23. Pre- 
vious work on the JCL model within VCA was carried 
out in Ref. 119|. 

A. Variational cluster approach for bosons 

The basic concept of VCA is that the grand potential 
Q is expressed as a functional of the self-energy Z and 
that Dyson's equation for the exact Green's function G is 
recovered at the stationary point of the self-energy func- 
tional rt[Y.]^^ In order to evaluate the self-energy 
Z of the investigated system is approximated by the self- 
energy of an exactly solvable, so-called reference, system. 
In practice, this means that the self-energy Z becomes a 
function of the set x of single-particle parameters of the 
reference system, i.e., Z = Z(x). For bosonic systems 
the approximated grand potential reads^i 

Q{x) = n\x) + Tr ln(-G'(x)) + Tr ln(-(Go ^ - Z(x))) , 

(4) 

where primed quantities correspond to the reference sys- 
tem and Go is the noninteracting Green's function. The 
stationary condition on Q(x) is given by 

OX 

This condition can be evaluated numerically by varying 
some or all of the single-particle parameters. In order 
to guarantee that a given physical quantity (such as the 
number of particles) is thermodynamically consistent, it 
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is necessary that Q is stationary with respect to the asso- 
ciated couphng constant (here the chemical potential). 
Therefore, varying uJc ensures that the total number of 
photons is thermodynamically consistent. On the other 
hand, it would be advisable for a conserved quantity, i. e., 
Np to be consistent as well. Otherwise uncommon situ- 
ations could occur. For example, as we show below, the 
total particle density Np/N ^ evaluated as a trace of the 
Green's functions is not integer in the Mott phase. This 
effect becomes stronger close to the tip of the Mott lobe, 
see Fig.[3](b). The noninteger particle density, occurring 
when fi is not taken as a variational parameter, clearly in- 
troduces an uncertainty in the determination of the phase 
boundary. 

In principle, however, there is a formal difficulty in 
taking /i as a variational parameter. The problem is re- 
lated to the coupling of /i with atomic excitations, which, 
in contrast to photons, cannot be seen as noninteracting 
particles. This is, in general, not allowed within VGA, 
whereby the reference system can differ from the physical 
one by a single-particle Hamiltonian only. The solution is 
readily overcome by observing that the two-level atomic 
system can be mapped onto a hard-core boson model. 
In this way, fi couples to the total number of "atomic" 
bosons plus photons, i.e., a noninteracting Hamiltonian. 
The hard-core constraint simply becomes a local (in prin- 
ciple infinite) interaction, which is common to the refer- 
ence and to the physical system. 

The mapping of the two-level excitations onto hard- 
core bosons is mathematically achieved by the following 
replacements 

Cr+ ^ h\ , (J- hi , 

i;,) ^ |0,) and It^) ^ . 

This is valid provided one excludes states with double 
occupation of h particles even as intermediate states.^^ 
With this mapping the JGL Hamiltonian reads 

^JCL = _ ^ ^ „J „ . _ A ^ bl b,+g ^(a, b] + a] h) 
-ill- We) Np + lim ^ ^'(^1 - 1) ' 

i 

(6) 

where we have formally implemented the hard-core con- 
straint by introducing an infinite interaction for b parti- 
cles. In the restricted Hilbert space of zero or one hard- 
core boson per lattice site, the matrix elements of the two 
representations are identical. In principle, states with 
higher occupation number blbi > 1 have to be consid- 
ered in the bosonic version as well. However, the occupa- 
tion of such states would cost infinite energy and, there- 
fore, they do not influence the energies obtained from the 
Hilbert space sector with occupation numbers blbi <li^ 
We have checked this aspect numerically for very large 
U. It can also be verified easily when the sector of the 
Hilbert space with blbi > 1 is included perturbatively. 



These considerations can be straightforwardly extended 
to light-matter models with more than one atom or atom- 
iclike structure (with two relevant levels) per cavity. In 
this case, one introduces a boson species for each atom 
and the hard-core constraint is enforced for each boson 
species. 

In our calculation we take both parameters cjc and e of 
the reference system as variational parameters (/i is just a 
linear combination), which ensures thermodynamic con- 
sistency for the particle number of both species, and, con- 
sequently, of the total particle number. We show below, 
that varying both parameters instead of just Uc provides 
an improvement in the accuracy of the phase boundaries 
for a given cluster size, see Tab. HI 

The present formulation of the VGA cannot address 
the superfluid phase, and is, thus, restricted to the Mott 
phaseJ^i2i Outside of the Mott lobes, peaks with neg- 
ative (positive) spectral weight appear in the positive 
(negative) uj region, signaling the instability towards the 
superfluid phase. A treatment of the superfluid phase 
requires a Nambu Green's function treatment analogous 
to the fermionic case. The boundary of the Mott phase 
could be, in principle, determined by this criterion. How- 
ever, it is simpler (and, of course, equivalent) to identify 
it, for a given t, as the region between the ground-state 
energies of the (A/^ + l) and (A/p — l)-particle states, which 
can be directly inferred from the single-particle spectral 
function. 

In VGA, the reference system is chosen to be a de- 
composition of the total system into identical clusters, 
which means that the total lattice of N sites is divided 
into clusters of size L. Mathematically this can be de- 
scribed by introducing a superlattice, such that the origi- 
nal lattice is recovered when a cluster is attached to each 
lattice site of the superlattice. The reference system de- 
fined on a cluster is solved by means of the band Lanc- 
zos method j^^i^^ The initial vector of the iterative band 
Lanczos method for the single-particle excitation term of 
the cluster Green's function contains 2L elements and is 
given by 

{a\ IV^o) , 4 \^o) ... 4 \^o) , crt \^o) . . . IV^o)} , (7) 

where Itpo) is the Np particle ground state. For the single- 
hole excitation term the initial vector of the band Lanczos 
method is obtained by replacing the creation operators 
in Eq. ([7]) by annihilation operators. 

To evaluate the grand potential and the single-particle 
Green's function of the original system we use the bosonic 
Q-matrix formalism. This formalism yields the Green's 
function G(k, uj) in a mixed representation, partly in real 
space and partly in reciprocal space, see App. [X] The 
matrix G(k, uj) is of size 2L x 2L and k belongs to the 
first Brillouin zone of the superlattice. Due to the specific 
order of the creation operators in the initial vector of the 
band Lanczos method we are able to extract the Green's 
function for photons G^^(k, u) and the Green's function 
for two- level excitations G^^(k, uj) from G(k, uj) in the 
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following way 

GP^,{k, a;) =G,,,(k, uj) and 

where r, s G [1 . . . L]. The application of the periodiza- 
tion prescription proposed in Ref. IsH (Green's function 
periodization) yields the fully k dependent Green's func- 
tions GP^{k, uj) and G^^(k, uj). From that we are able 
to evaluate the single-particle spectral function 



-_ImG^(k, uj) 

TT 



A^k, uj) 
the density of states 

N^iuj)^ j A^{k, a;)dk=l^A-(k, uj) 

and the momentum distribution 



(8) 



(9) 



n^(k) = - j A^(k, uj)duj , 

•^^ — OO 



(10) 



where x can be either ph for photons or ex for two-level 
excitations. We use the Q-matrix formalism to evaluate 
the momentum distribution, since this approach yields 
particularly accurate result s.^^ Furthermore we calculate 
the spatial correlation functions 



C^il' = {A^j) and 



(11) 



which just depend on the distance between two cavities i 
and j, i.e., Cf^ = C^(|r^ — r^l). Notice that the poles of 
the hard-core boson Green's function coincide with the 
poles of the two-level excitation Green's function as the 
energies of both representations are identical. However, 
the hard-core boson Green's function exhibits additional 
poles located at energies of the order U ^ oo. The addi- 
tional poles which have finite weight result from the fact 
that excitations such as h\ are in principal allowed 
but cost infinite energy, whereas the corresponding exci- 
tation (j^ Wi) is strictly forbidden. Therefore, the single- 
particle correlation functions {h\^{t)h\^ and {cF^{t)cF^) 
differ only by contributions from frequencies of the order 
U ^ (X). Yet it should be mentioned that the single- hole 
correlations function of hard-core bosons is not affected 
by these considerations as {h\it) h\^) is always equiva- 
lent to (cTj^ {t) cTj^ ) . This also implies that the spectral 
weight of the poles with negative energy are identical for 
both representations and that the particle density of the 
two-level system is equal to the particle density of the 
hard-core bosons. In the following, we will always speak 
loosely about two-level excitation Green's functions but 
we have to keep in mind that there are differences in the 
single-particle spectral weight of the hard-core boson and 
two-level excitation Green's functions at infinite energies. 



B. Polariton properties of the quasiparticles 

In the next step, we want to investigate the polar itonic 
properties of the JCL model, which arise due to the cou- 
pling between the photons and the two- level excitations. 

Adding a particle or hole to the many-body ground 
state may result in quasiparticle or collective excitations 
which are built up by the {Np ± l)-particle eigenstates 
of the many-body system entering the Green's function. 
These many-body eigenstates for the infinite system can 
be extracted within the VGA framework from the VGA 
Green's function. As shown in App. [Aj they are linear 
combinations of the particle and hole excitations of the 
cluster Green's function weighted by the eigenvector ma- 
trix X, defined in App. [XI 

Our goal is to describe the eigenvectors of the (A/p±l)- 
particle Hilbert space, which form the quasiparticle exci- 
tations of the Green's function by polaritonic quasiparti- 
cles added to the exact Np particle groundstate IV^o). To 
this end, we introduce the polariton creation operators 
^1 k particle excitations and h^^ ^ for hole excitations 
as appropriate linear combinations of photons and two- 
level excitations 

Pi,k = /3p(k)4 + 7;(k)a+ , (12a) 
hl^ = Pt{^^)a^ + ^t{\^)a^ . (12b) 

It should be stressed that the hole creation operator 
is not the adjoint of the particle creation operator or its 
annihilation counterpart, which it would be in the case of 
noninter acting particles. As we will see, the coefficients 
or weights of the linear combinations P^j^ik) and 7^^^(k) 
depend on the wave vector k, the quasiparticle band a, 
and additionally on the filling n, which is not explicitly 
written in Eq. ([12]), since the filling dependence is not 
important for the present discussions. The normalized 
polariton quasiparticle states are defined by applying the 
polaritonic operators on the exact Np particle ground 
state l^/^o) yielding 



Km) 



Y''{'*/'o|/la,k/ll,kl^o) 



and 



(13a) 
(13b) 



respectively. The normalization terms can be rewritten 
as 

(^oba,kpl,kl^o) = z^^(k)Sp(k)z^(k) and (14a) 
(V'ol Km l^Po) = z^^(k) S,,(k) z^(k) . (14b) 

In Eq. ([T4|) the vectors z^^^(k) are defined as z^y^(k) = 
(/3^/^(k), 7^/^(k))^ and Sp/^(k) are the overlap matrices 
of single-particle excitations and single-hole excitations. 
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respectively. The overlap matrix for the hole excitations 
is given by 



S/.(k) 



(4«k) (4^k) 



where the static correlation functions are evaluated in the 
Np particle ground state IV^o)- All quantities entering Sh 
are correctly evaluated in the hard-core boson model as 
no excitations of the "two-level bosons" into the n > 1 
sector occur. For the particle case the situation is differ- 
ent, as we need to evaluate 



Sp(k) 



(^k «l) 



The term (cr^cr^) of the two- level system cannot be di- 
rectly evaluated in the hard-core boson model. Using the 
commutator property [cr^~,cr^] = for i 7^ j and the lo- 
cal anticommutation relation {cr~^cr^} = 1, we end up 
with an expression that only contains static correlation 
functions which can be computed correctly within the 
hard-core boson model 



S,(k) = ( 



l + (^^^k)-|rEkK<^k: 



In order to derive a formalism to construct the opti- 
mal polariton weights, we start out with the analysis of 
an exact eigenvector IV^^^^^) of the Hamiltonian in the 
{Np + l)-particle sector. For the sake of clarity we will 
suppress in the following considerations the index k for 
all quantities, and the indices a and p for quasiparti- 
cle weights and wave functions. The optimality criterion 
in this case is clearly the overlap of the exact eigenvec- 
tor with the approximate (normalized) vector given in 
Eqs. ([ni) and (p]) 



where / denote the components of the two-dimensional 
vectors, and di^k = <^k and (i2,k = cr^ , see App. [A] The 
maximization of | (V^^^^ l^/^j,) p leads to the generalized 
eigenvalue problem 

A"z" = ASpz" , (15) 
where the elements of the 2x2 matrix A^ are 



In Eq. ([T5j) we replaced z^ by z^ as the eigenvalues are 
just determined except for a constant Z, which will be 
specified later. As the eigenvalue corresponds to the 
value of the overlap squared A = | (ipj^'^'^lipiy) p, the de- 
viation of the eigenvalue from one is a measure of the 
quality of the polariton approximation. It also points out 



that the eigenvector corresponding to the largest eigen- 
value determines the optimal polariton coefficients. In- 
terestingly, A^j is the contribution of the excitation u to 
the corresponding spectral function, i. e., its quasiparticle 
weight. In general, the quasiparticle peak is a superposi- 
tion of several exact many-body eigenstates. Hence, the 
obvious generalization of the optimality criterion is to 
sum over all eigenstates which contribute to the quasi- 
particle excitation a. To this end we define an energy 
window Qa in which the quasiparticle peak a is located 
and we integrate the spectral density in this energy win- 
dow resulting in 



i/j(k, Qa) 



The polariton coefficients are again obtained by the gen- 
eralized eigenvalue problem 



A(k, Qa)^ = 
and the eigenvalue is given by 



A Sp z 



A 



zf A(k, 
zt So z 



(16) 



The eigenvalues are still restricted to the unit interval 
[0, 1]. The lower limit is due to the positivity of A and Sp. 
The upper limit follows from the property that a summa- 
tion of the integrated spectral density over all nonover- 
lapping energy intervals is given by 

particles 

E ^j(k, fta) = {di,\.dl^^) = {Sp)ij . 



Of course, z and, hence, the polariton operators will de- 
pend on the wave vector k, the quasiparticle band index 
a and the filling n, i. e., the Mott lobe. The discussion so 
far was for the particle case only, however, it is straight- 
forward to iterate the procedure for the hole case. 

Eventually, we merely need the integrated spectral 
density A(k, Q^) determined within the VGA framework, 
which is given by 

Aij{k,na) = - Yl (QX)/,.(x-iSQt),,j . 

Details are presented in App. [A] as well as the proof that 
all contributions of the sum have the same sign, which is 
necessary for the optimality criterion to make sense at all. 
The optimality criterion as well as the eigenvalue problem 
only fix the coefficient vector z up to a normalization 
factor Z, i.e., z = Z z. The latter is determined by 
the condition that the total spectral weight should be 
conserved 



Z^z^Kz = tr A . 



(17) 



As the excitations can now be described by wave vec- 
tor, band and filling dependent polaritonic quasiparti- 
cles, it remains to evaluate the polariton spectral function 



6 




FIG. 2: (Color online) Phase boundaries of the JCL model 
in one dimension for zero detuning A = 0. (a) VGA results 
for the variational parameters x = {cjc, e} and various cluster 
sizes of the reference system. The gray shaded area indicates 
DMRG data."^ (b) Phase boundaries obtained for the largest 
cluster (L = 8 for the first Mott lobe and L = 6 for the second 
Mott lobe). The marks refer to parameters where spectral 
functions are evaluated. 



A^(k, cj), which is due to the invariance of the trace in 
Eq. (fT7|) equal to the sum of the photon spectral function 
A^^(k, (ju) and the two- level excitation spectral function 



IV. RESULTS 

In this section, we present the results of our calcula- 
tions. Specifically, in Sec. IIV A| we discuss the quantum 
phase transition from Mott phase to superfluid phase oc- 
curring in the JCL model and investigate the impact of 
the variational parameter space on the accuracy of the 
results. In Sec. IIVB| we study the spectral properties of 
both photons as well as two- level excitations. The first 
two subsections refer to results obtained for zero detun- 
ing A = 0, whereas nonzero detuning is considered in 
the third subsection. Sec. IIV CI Finally, in Sec. IIVD| we 
study the polaritonic properties of the JCL model. In 
particular, we introduce polariton quasiparticles as wave 
vector and filling dependent linear combinations of pho- 
tons and two-level excitations and analyze the weights of 
their constituents. 



A. Quantum phase transition 

The JCL model exhibits, comparable to the BH 
model, ^ a quantum phase transition from a localized 
Mott phase to a delocalized superfluid phase. For in- 
teger particle density and small hopping strength t, the 
ground state of the system is a Mott state. The first two 
Mott lobes of the one-dimensional (ID) JCL model for 
zero detuning A = obtained by means of VC A with the 
variational parameters x = {cjc, e} are shown in Fig. [21 
As discussed in the previous section, including e in the set 
of variational parameters is nontrivial and is solely pos- 
sible since the two-level excitations can be mapped onto 




t t 

(a) (b) 

FIG. 3: (Color online) Comparison between the results ob- 
tained with the variational parameters x = {cjc, e} and 
X = {oJc}, respectively, for small clusters of size L = 2 and 
L = 4. (a) Phase boundaries at the tip of the first Mott lobe. 
The gray shaded area indicates DMRG results.— (b) Total 
particle density n, which is the sum of the photon density and 
the two-level excitation density, across the first Mott lobe. 

hard-core bosons. The gray shaded area in Fig.[2](a) in- 
dicates DMRG results for the phase boundary obtained 
by D. Rossini et al in Ref. |T3. We find excellent agree- 
ment between the phase boundary evaluated by means 
of VCA with the variational parameter set x = {cjc, e} 
and the DMRG results, even at the lobe tip, where quan- 
tum fluctuation effects are most important, and even for 
moderate cluster sizes L ^ 4. Figure [3](a) compares the 
phase boundaries at the tip of the first Mott lobe for dif- 
ferent variational parameters. The results obtained with 
X = {cjc, e} are connected by lines, whereas the open 
symbols correspond to x = {cjc}. We observe that us- 
ing both on-site energies as variational parameters im- 
proves the results for the phase boundary and also yields 
a better approximation for the slope of the lobe tip. A 
quantitative measure for the quality x of the calculated 
phase boundary is given by the absolute deviation from 
the DMRG data per phase boundary point 

x = l4;EW'-P?|, (18) 

where pf and are corresponding phase boundary 
points calculated by means of VCA and DMRG, respec- 
tively, and Mp is the number of phase boundary points, 
which contribute to the sum. In Tab. H] we compare the 
quality x/10~^ of the phase boundary between the two 
sets of variational parameters for various cluster sizes. 
When using the augmented set of variational parame- 
ters X = {cjc, e} in contrast to x = {oJc} we observe an 
improvement in the quality of the phase boundary which 
ranges from 1.3 to 1.7 depending on the cluster size of the 
reference system. Using both the resonance frequency ujc 
of the cavities and the energy spacing e of the two-level 
system as variational parameters thus provides a signif- 
icant improvement with respect to the case of a single 
variational parameter.^- As discussed in Refs. |27| and[30, 
a correct particle density in the original system can only 
be obtained when the corresponding on-site energies are 
included in the set of variational parameters, i.e., in the 
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TABLE I: Quality x/10 ^ of the phase boundary for x = {ujc} 
and X — {(jjc, e}, respectively. The quality x is evaluated using 
Eq.dEl). 

L . . . number of cluster sites 
variational parameters 

improvement in quality when using the variation- 
al parameters x = {cjc, e} instead of x = {cJc} 



IMP 



L 




Wc, e} 


IMP 


2 


15.95 


11.34 


1.41 


4 


8.20 


4.92 


1.67 


6 


5.34 


3.16 


1.69 


8 


3.95 


3.07 


1.29 



case of the JCL model x = {cjc, e}. This is demonstrated 
in Fig.[3](b), where the total particle density n, which 
consists of a photon and a two-level excitation contribu- 
tion, is evaluated along the first Mott lobe. For x = {oJc} 
the deviation of the particle density from one is growing 
with increasing hopping strength t but shrinking with in- 
creasing cluster size L. However, when e is included as 
variational parameter the total particle density n is as 
desired equal to one across the whole first Mott lobe. A 
deviation of about 0.001 can be observed for t = 0.2. Yet, 
the hopping strength t = 0.2 is probably even slightly 
above the critical hopping strength t*, which indicates 
the tip of the Mott lobe 

The phase diagram of the ID JCL model is in many 
aspects similar to the phase diagram of the ID BH 
model i^^i^^ Particularly, the Mott lobes are point shaped 
and a reentrance behavior can be observed, which means 
that for certain values of /i upon increasing t the system 
leaves the Mott phase and later on enters it again. Yet a 
very important difference is that the width of the lobes of 
the JCL model at zero hopping is shrinking with increas- 
ing particle density. This comes from the fact that the 
effective on-site repulsion of the JCL lattice model, which 
is hidden in the interaction between photons and two- 
level excitations, is not constant, as in the Bose-Hubbard 
model. The exact location of the phase boundaries at 
zero hopping is derived as a by-product in App.[Bl whose 
major intention is, however, to introduce the notation 
used for the dressed states |n, a) and for the correspond- 
ing energies E^^^a)^ where a G { — , +} describing the 
ground state and the excited state in the correspond- 
ing constant particle number sector of the single-cavity 
Hilbert space. 
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FIG. 4: (Color online) Photon spectral function A^^(k, cj), 
first row, and density of states A^^^(cj), second row. Two- 
level excitation spectral function A^^(k, cj), third row, and 
density of states N^^(uj), fourth row. The spectral functions 
are evaluated for the parameters (a) t = 0.03, ii — ujc = —0.75, 
A = and (b) t 0.12, - cjc = -0.84, A = 0, which be- 
long to the first Mott lobe. The dashed lines in the spectral 
functions in (a) correspond to first-order degenerate pertur- 
bation theory results, see App. [C] The Roman numerals in 
the captions of the subfigures refer to the marks in Fig.[2](b). 



B. Spectral properties of photons and two-level 
excitations 

The spectral function for photons A^^(k, a;), the spec- 
tral function for two- level excitations A^^{k^ uj) and the 
corresponding densities of states N'^^iuo) and N^^iuj) 
evaluated by means of VCA for parameters belonging 
to the first Mott lobe are shown in Fig.HJ We use an 



artificial broadening r] = 0.03 and the variational pa- 
rameter set X = {cjc, e} for the numerical evaluation of 
the spectral functions. Both spectral functions 74^^(k, uo) 
and 74^^(k, uo) have the same gap as the photons and the 
two-level excitations are coupled. The spectral functions 
of the JCL model generally consist of four bands. This 
can best be understood in terms of the analytic solution 
of the JCL model for zero hopping strength t = 0. The 
ground state IV^o) of the JCL model in the Mott phase 
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with particle density n for zero hopping is given by the 
tensor product state 

N 

l*o) = (g) |n, , (19) 

iy=l 

where |n, — )^ is the dressed n particle ground state of 
lattice site u. The states with a single-particle excitation 
are those, where N — 1 sites remain in the dressed state 
|n, — ) and one site is excited to the state \n + 1, a). Sim- 
ilarly, for the single-hole excitation N — 1 sites remain 
in the state |n, — ) and one site is excited to the state 
|n — 1, a). In both cases, the excited states are TV fold 
degenerate as the particle/hole excitation can be located 
on any of the N lattice sites. The degenerate states have 
thus the structure 

N 

1^^' ^) = |n + 1, a), (g) |n, and (20a) 

1^=1 

|*«'')^|n-l,a>,(g)K-), , (20b) 

u = l 

respectively. Two of the four bands, we refer to them 
as lower modes emerge from the excitation of site 

i from the dressed state |n, —)■ to the states |n ± 1, — )^, 
which are ground states of the corresponding Hilbert- 
space sector with constant particle number. Analogously, 
we refer to the bands which emerge from the excitation 
of site i from |n, —)■ to the excited states in the corre- 
sponding particle sector |n ± 1, as upper modes oo^^j^- 
The presence of the upper modes has been first noted by 
S. Schmidt et al in Ref. |16| and has been numerically 
observed in latest QMC calculations^^ as well. The two 
upper modes (^^^^ indicate a clear deviation from the BH 
physics, which emerges due to the composition of two dis- 
tinct particles. As discussed in the previous section, the 
two particle bands o;^, a G { — , +}, determine the po- 

lariton particle creation operators ^ whereas the two 

hole bands uj^ specify the hole creation operators /ij^ ^. 

In the spectral functions of Fig.[4l the lower modes 
^p/h correspond to the cosinelike shaped bands centered 
around = 0. The intensities of the lower modes o;"^^ 

are contrary for the photon spectral function 74^^(k, uj) 
and the two- level excitation spectral function A^^Ck^uj). 
For 74^^(k, u) the particle band uj~ is more intense than 
the hole band uj^ whereas the hole band is more intense 
than the particle band for A^^(k, cj). For the first Mott 
lobe the upper hole mode u'^ does not exist as this would 
require to excite a single-site i from the dressed state 
|1, — )• to the non-existing state |0, +)-. Thus, only the 
upper particle mode cj+ can be observed in the spectral 
functions shown in Fig. 21 which corresponds to the es- 
sentially flat band located at — /i ~ 3. In App. [Cl 
we evaluate the single-particle and single-hole excitation 
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FIG. 5: (Color online) Extract of the lower particle band cj" 
of the photon spectral functions A^^(k, u) shown in Fig.[4l 
where the parameters (a) t — 0.03, ii — ujc = —0.75, A = and 
(b) t = 0.12, /i-cjc = -0.84, and A = have been used. VGA 
results (density plot) are compared with bands evaluated by 
means of first-order degenerate perturbation theory (dashed 
lines) . 



bands by means of first-order degenerate perturbation 
theory, which yields 

ujp^i = {ujc — /i) + a q{n + 1) + q{n) — 21^ cos k and 

(21a) 

1 = (^c — iJi) — OL q{n — 1) — q{n) -\-2t^ cos k , (21b) 

respectively, where t^^^ is the renormalized hopping 
strength. Figure [4](a) shows, additionally to the spectral 
functions obtained by means of VGA, the perturbation 
results for the bands. For small hopping strength we 
observe, as expected, good agreement between the two 
approaches. From the analytic solution of the bands we 
are able to extract their width, which is given by 2t^^^. 

The renormalization factor in i^^^ essentially consists of 

a square of the form (a + 6)^, see Eqs. (|C5p and (|C7p . 
Evaluating these expressions shows that a, 6 > for the 
lower modes but a > and 6 < for the upper 

modes ou^^^- Therefore, a and b almost cancel each other 
in the latter case, which yields a small renormalized hop- 
ping strength of the upper modes i^^j^ in comparison to 

the one of the lower modes and thus, essentially flat 
upper particle/hole bands uj^^^t^ Plugging in the value 
of the modified chemical potential ii — ujc = —0.75, which 
has been used to evaluate the spectral function shown 
in Fig.[l(a), into Eg. (f2Ta|) yields cj^^ ^ 3.16, where 
we neglected the dependence on the wave vector. This 
matches perfectly with the VGA results. In addition to 
previous worki^*^ we evaluate the upper modes not only 
for photons but also for two level-excitations. Interest- 
ingly, the spectral weight differs significantly for the two 
types of particles. In particular, the upper particle mode 
uj^ has a very large intensity in the two-level excitation 
spectral function 74^^(k, cj), but is almost not visible in 
the photon spectral function A^^{\^^ uj). For the spectral 
function shown in Fig. [41(b) a different chemical poten- 
tial fi — CjOc = —0.84 has been used. Thus, the upper 
particle mode is shifted slightly upwards in comparison 
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— VGA 
— cluster 




FIG. 6: (Color online) Comparison between the density of 
states obtained from the VCA Green's function, solid lines, 
and the density of states obtained from the cluster Green's 
function, dashed lines, (a) density of states of photons N^^{u) 
and (b) density of states of two- level excitations N^^ {uj). The 
parameters used for these plots are the same as in Fig.[4](b). 



to Fig.|4](a) and is located at ( 



3.25. Figure [5] shows 



the lower particle band uo^ of the photon spectral func- 
tion 74^^(k, ijj) for the same parameters as in Fig.lH In 
this figure we compare the VCA results for different hop- 
ping strengths with the results obtained by means of first- 
order degenerate perturbation theory. For smah hopping 
strength, t = 0.03, see Fig.[5](a), the perturbative results 
agree very well with the VCA results in both the width 
as well as the shape of the band. However, for large hop- 
ping strength t = 0.12, which is already close to the tip 
of the Mott lobe, the lower particle band does not ex- 
hibit a simple cosine shape anymore, see Fig.[5](b). In 
addition the width of the band is slightly overestimated 
by first-order degenerate perturbation theory. 

In the spectral functions shown in Fig.|4](b) there is 
additional spectral weight located at — /i ~ 2. We can 
exclude that this additional weight stems from the pe- 
riodization prescription used in VCA or from any other 
VCA internal processes as it also appears in the cluster 
Green's function, which is solved by exact diagonaliza- 
tion. This can be verified best by comparing the density 
of states obtained from the VCA Green's function with 
the density of states obtained from the cluster Green's 
function, see Fig.O Both densities of states, the one 
obtained from the cluster Green's function and the one 
obtained from the VCA Green's function, exhibit a peak 
located at — /i ^ 2. The additional peak can be re- 
vealed in the framework of perturbation theory. First- 
order local particle fluctuations in the ground state will 
have contributions of the form 



N 



(1) 



(g) n ■ 



In, 



1^ = 1 



where V correspond to nearest-neighbor sites. Due to 
the energy denominator the predominant terms are 
those with a = /3 = — . The correction term |A?/;^^^) is 
proportional to the hopping strength t, which explains, 
why the additional peak is not present in Fig.|3](a). The 
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FIG. 7: (Color online) Photon spectral function A^^(k, cj), 
first row, and density of states N^^{(jj), second row. Two- level 
excitation spectral function A^^(k, cj), third row, and density 
of states N^^(u), fourth row. The spectral functions are eval- 
uated for the parameters (a) t = 0.002, ii — ujc = —0.37, A = 
and (b) t = 0.012, fi - Uc = -0.38, and A = 0, which be- 
long to the second Mott lobe. The Roman numerals in the 
captions of the subfigures refer to the marks in Fig.[2](b). 



particle excitation couples to final states with an addi- 
tional particle either on site l^V oi on one of the remain- 
ing sites. A detailed analysis shows that the excitation, 
responsible for the additional peak at about a; — /i ^ 2, 



N 

|V-^.+i) = |n+l, -),®\n, +),, (g) K . 
The corresponding excitation energy is given by 
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FIG. 8: (Color online) Momentum distribution (a) in the first 
Mott lobe and (b) in the second Mott lobe for the parameters 
marked with Roman numerals in Fig.[2](b). Solid lines cor- 
respond to the momentum distributions of photons n^^(k) 
and dashed lines to the momentum distributions of two-level 
excitations n^'^(k). 

= cjc - M - q{n + 1) + ?>q{n) . 

For zero detuning and ji — uOc = —0.84 the energy is 
Cbp = uo — II = 2.4. 

As discussed before the upper hole mode uo^ does not 
exist in the first Mott lobe. Yet, the mode uj^ is present 
in spectral functions of the second Mott lobe, see Fig. [71 
According to Eq. (|2T]) the upper modes are located at 
^pi ^ ^-^^ '^hi ^ —2.04 for the parameters used 
in Fig. [71(a). This matches very well the results obtained 
by means of VGA. The chemical potential of the spectral 
function shown in Fig. [71(b) differs from the one of (a) 
merely about 0.01. Thus, the bands cj^//^ are located at 
rather the same position in both spectral functions. 

The momentum distribution for photons n^^(k) and 
two- level excitations n^^ (k) in the first and second Mott 
lobe are shown in Fig. [51 For increasing hopping strength 
t the momentum distribution becomes more peaked for 
both the photons and the two- level excitations. In the 
first Mott lobe the momentum distributions n^^(k) and 
n^^(k) are centered around 0.5, which means that the 
cavities are on average equally occupied by photons and 
two- level excitations. In the second Mott lobe n^^(k) is 
centered around 1.5. However, n^^(k) is still centered 
around 0.5, as the maximum local occupation number of 
the two-level systems is restricted to one. 

In order to display the slowing down of correlations 
upon approaching the boundary of the Mott phase, we 
evaluate the spatial correlation function C^(|r^ — r^l) in 
the first Mott lobe (Fig.[9|). The spatial correlation func- 
tion can be obtained from the Fourier transform of the 
momentum distribution. For small distances Ir^ — r^l 
between sites i and j the correlation function is a super- 
position of multiple exponential functions with distinct 
strengths of decay. For large distances, however, the ex- 
ponential function with the smallest decay dominates and 
thus the correlation function is of the form 

C^(|r,-r,-|)oce-""l"----^l , (22) 



FIG. 9: (Golor online) Correlation function (a) for photons 
and (b) for two-level excitations in the first Mott lobe. The 
Roman numerals in the legend refer to the parameters marked 
in Fig.[2l(b). 
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FIG. 10: (Golor online) Phase boundaries of the ID JGL 
model for the detuning (a) A = — 1 and (b) A = 1. The 
marks refer to the parameters where spectral functions are 
evaluated. 



as expected in the insulating phase. Using VGA we are 
able to extract the correlation length = 1 /a^ , as data 
are available for large distances between two sites i and 
j. From a linear fit for sufficiently large distances we 
obtain = a^^ = 1.711 ± 0.001 for the parameters I, 

see marks in Fig.[2l(b), and = aff = 0.317 ± 0.001 
for the parameters II. Therefore, the slope of the corre- 
lation function is the same for the two particle species, 
which is due to the coupling between the photons and 
the two-level excitations. As in the BH model^^ the ab- 
solute slope of the correlation function shrinks with 
increasing hopping strength, which is a precursor of the 
superfluid phase, where the correlation between sites per- 
sists up to long distances. 



C. Nonzero detuning 

The detuning A, which is the difference between the 
resonance frequency uOc of the cavities and the energy 
spacing e of the two- level systems, is a very important 
parameter of the JGL model. By varying the detun- 
ing it is possible to change the width of the Mott lobes. 
Phase boundaries obtained by means of VGA with the 
set of variational parameters x = {ujc^ e} for A = — 1 
and A = 1 are shown in Fig.[TOl For the parameters 
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FIG. 11: (Color online) Photon spectral function A^^(k, cj), 
first row, and density of states N^^{(jo), second row. Two- 
level excitation spectral function A^^(k, cj), third row, and 
density of states N^^(lu), fourth row. The spectral functions 
are evaluated for the parameters (a) t = 0.036, — ujc = 
-0.54, A = -1 and (b) t = 0.2, /i - cj^ = -1.2, A = 1. The 
Roman numerals in the captions of the subfigures refer to the 
marks in Fig. 1101 



marked with x we evaluate the spectral function of pho- 
tons 74^^(k, uj) and two- level excitations 74^^(k, a;), see 
Fig.[TTJ An interesting effect can be observed in the 
spectral functions 74^^(k, uo). Namely, the intensity of 
the upper band a;+ depends significantly on the detun- 
ing A. For negative detuning A = — 1, the upper mode 
in 74^^(k, uj) is very intense, see Fig.[TT](a), whereas it 
is almost not visible for positive detuning A = 1. This 
behavior remains valid when the spectral functions for 
positive and negative detuning are evaluated for identi- 
cal hopping strength. The zero-hopping result for the 
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FIG. 12: (Color online) Momentum distribution evaluated for 
the parameters marked in Fig. llOl where (a) corresponds to 
the parameters V, i. e., negative detuning A = — 1 and (b) to 
the parameters VI, i.e., positive detuning A = 1. 



energy of the upper mode is uj^^^ ~ 3.15 for the spec- 
tral function shown in Fig.[TT](a) and uo^^^ ^ 3.82 for the 
spectral function shown in Fig.[TT](b). The momentum 
distributions of photons n^^ (k) and two- level excitations 
n^^(k) for the parameters marked in Fig. [10] are shown 
in Fig. [121 For negative detuning it is energetically more 
expensive to excite the two-level system than to add a 
photon to the cavity. Thus, the momentum distribution 
of photons n^^(k) dominates over the momentum distri- 
bution of two- level excitations n^^(k). For positive de- 
tuning the situation is reversed and n^^ (k) is larger than 
n^^(k) for all values of the momentum. 



D. Polariton quasiparticles 

Up to now we investigated the photon properties and 
the two-level excitation properties of the JCL model sep- 
arately, by extracting the Green's function of photons 
G^^(k, uj) = G^^^^ {uj) and the Green's function of two- 
level excitations G^^(k, uo) = G - +{uj) from the com- 

k k 

pound Green's function G(k, ui), which is a 2 x 2 matrix 
of the form 



G(k, w) = 



G t(w) G 



G - t(w) G-^+(lo) 



(23) 



Next we will discuss the polaritonic properties of the 
JCL model. We start out with the first Mott lobe for zero 
detuning and focus again on the parameter set marked 
as II in Fig.[2l i. e., t = 0.12, /i - cj^ = -0.84 and A = 0. 
The polaritonic spectral function A^(k, cj) and the cor- 
responding density of states N^{uj\ which is by con- 
struction identical to the total density of states of pho- 
tons plus two-level excitations, is shown in Fig. [131 For 
the first Mott lobe the hole case is special since both, 
G~ |n, — ) oc |0, — ) and a |n, — ) oc |0, — ) yield the exact 
zero-particle state. Consequently, the polariton can be 
chosen ad libitum, it will always be exact. Therefore in 
Fig. [14] only the particle part of the polaritonic weights 
is depicted. The right panel represents the result for the 
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FIG. 13: (Color online) Polariton spectral function (a) and 
density of states (b) evaluated for the parameters II corre- 
sponding to the first Mott lobe, i.e., t = 0.12, /a — ujc = 
—0.84, and A = 0. In (b) the polariton density of states 
N^{(jj) is compared with the sum of the photon density of 
states A^^^(k, u) and the two- level excitation density of states 
A^^^(k, cj), which coincide by definition. 



! 0.5 



(a) 



— P — y ---X 



0.5 



(b) 



-0.5 





k/7i 



0.5 



-0.5 





k/7i 



0.5 



FIG. 14: (Color online) Photon contribution [3 and two-level 
excitation contribution 7 to the polariton quasiparticles. (a) 
shows the results for pj^ corresponding to the upper particle 
band uOp and (b) for pl_ ^ corresponding the lower particle 
band cj" . Additionally to the weights /3 and 7 the overlap A 
is shown. 



lower particle excitation. The polariton has very pro- 
nounced photonic character and the weights of photons 
and two-level system have opposite sign. Interestingly, 
the lower particle excitation can very well be mimicked 
by a single polariton on top of the A^^-particle ground 
state, as can be inferred from the fact that A ^ 1. More- 
over, a slight k-dependence of the weights is observed. 
Contrarily in the upper particle band, the polariton has 
pronounced two- level-system character, the weights have 
the same sign, there is almost no k-dependence, and the 
polariton description is poor (A ~ 0.2). 

Now we turn to the second Mott lobe, which allows 
us to study the hole polariton as well. The polariton 
spectral function and the corresponding density of states 
evaluated for the parameters IV, i. e., t = 0.012, ii — uOc = 
—0.38, and A = 0, are shown in Fig. [151 The weights 
are shown along with the overlap A in Fig. [161 The lower 
bands cj"^^ are well described by the quasiparticles as 
the overlap A is almost one for both bands. The upper 
bands co'^/^, however, are not described that well. In 
particular A ~ 0.2 for the upper particle band and A ~ 
0.85 for the upper hole band. The weights {3 and 7 are 



FIG. 15: (Color online) Polariton spectral function (a) and 
density of states (b) evaluated for the parameters IV corre- 
sponding to the second Mott lobe, i.e., t = 0.012, /a — ojc = 
-0.38, and A = 0. 



significantly more wave vector dependent, especially for 



the upper bands 



p/h 



, 1. e., a 



-. Apart from the more 
pronounced k-dependence, the weights for the particle 
case are rather similar to those of the first Mott lobe. 
However, there are striking differences in the weights for 
the particle and hole part within the second Mott lobe. 
First, the k dependence is more pronounced. Second, the 
sign of the relative weights is positive for both bands a = 
±, and finally, the composition of the polariton in the 
two bands is reverse. The lower band has predominantly 
photonic character, while opposite holds for the upper 
band. 

Eventually, we want to compare the VGA results 
with those of the single-site problem, which are de- 
rived in App. [D1 In the single-site problem the sign 
of the relative polaritonic weights is the same as that 
observed in the lattice. In the first Mott lobe the rela- 
tive weights for the particle case are for the upper band 
= 7^n=i//^^n=i = V2 + 1 and for the lower band 
the reciprocal relation holds q- = Pp^n=i/^p,n=i ~ ~Q+- 
There is agreement in the relative signs and the compo- 
sition of the polariton between the single-site limit and 
the lattice system, but the reciprocal property is strongly 
violated in the lattice case. This might be understood as 
follows. The itinerant particles are the photons. In order 
to gain kinetic energy it is convenient for the system to 
increase the photonic character in the dispersive lower 
band, depicted in Fig.[T5l(a). The upper band, on the 
other hand, has little dispersion and behaves more like 
the single-site limit. 

In the second Mott lobe, the relative weights obtained 
in the single-site limit for particle excitations are (7+ = 
VS -\- V2 and q- = —q-\-. Like in the lattice case, the 
weights of the particle part are comparable in the first 
and second Mott lobe. Quantitatively, the relative weight 
\q± \ is roughly 30% larger in the second Mott lobe, which 
is also the case in the lattice system. As far as the hole 
part is concerned, the single-site limit nicely corroborates 
all observations of the lattice model. 

In the single-site problem, the exact many-body eigen- 
states |n ± 1, a) can be generated exactly by suitable po- 
lariton operators acting on the state |n, — ). This is no 
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FIG. 16: (Color online) Photon contribution and two-level 
excitation contribution 7 to the polariton quasiparticles. (a) 
shows the results for corresponding to the upper particle 
band cj^, (b) for pl_ ^ corresponding the lower particle band 
Up , (c) for hl_ ^ corresponding to the lower hole band cj^ 

and (d) for /ij^ ^ corresponding to the upper hole band cj^. 
Additionally to the weights /3 and 7 the overlap A is shown. 



longer the case in the lattice due to local particle num- 
ber fluctuations induced by particle motion. Already in 
the single-site limit, the polariton operators are, how- 
ever, not universal, they depend on the filling n and in 
the lattice case even on the wave vector k. On top of 
that, the polariton operator for holes is not the adjoint 
of the corresponding polariton creation operator of the 
particle type, or in other word its annihilation operator. 



V. CONCLUSIONS 

In this paper we presented and discussed the spec- 
tral properties of the Jaynes-Cummings lattice model in 
one dimension obtained within the variational cluster ap- 
proach. Using the resonance frequency Uc of the cavities 
and the energy spacing e of the two-level systems as vari- 
ational parameters in the variational cluster approach 
procedure provides a significant improvement with re- 
spect to the case of a single variational parameter. On 
the one hand, varying both Uc as well as e (or, at least fi) 
is necessary to achieve a correct particle density in the 
original system and on the other hand improved results 
for the phase boundaries, and thus, for the spectral func- 
tions as well, are obtained due to the augmented set of 
variational parameters. In order to apply the variational 
cluster approach and include e as variational parameter 
the two-level systems have been mapped onto hard-core 



bosons, which yields correct poles of the Green's func- 
tion in the relevant energy range. We evaluated and dis- 
cussed spectral functions for photons and two-level exci- 
tations. The spectral functions generally consist of four 
bands, cosinelike shaped lower particle/hole bands, which 
are centered around zero energy, and essentially flat up- 
per particle/hole bands. An exception are the spectral 
functions in the first Mott lobe, which contain the two 
lower bands but from the upper bands only the parti- 
cle part. Using first-order degenerate perturbation the- 
ory, we evaluated analytical expressions for the bands, 
which allowed us to explain why the upper modes are 
essentially fiat whereas the lower modes exhibit a pro- 
nounced cosinelike shape. Additionally, we compared 
the analytical solution for the bands with the variational 
cluster approach results. For small hopping strength t 
we observe, as expected, good agreement between the 
two approaches. However for parameters located close 
to the tip of the Mott lobe, first-order degenerate per- 
turbation theory yields results that differ from the exact 
ones in both, shape and width of the bands. Further- 
more, we evaluated densities of states, momentum dis- 
tributions and spatial correlation functions for photons 
and two- level excitations. We also investigated detun- 
ing effects on the spectral properties and found indica- 
tions that the intensity of the upper particle band of the 
two-level excitation spectral function depends strongly 
on the detuning. Based on the information obtained 
from the photons and two-level excitations we investi- 
gated the polaritonic properties of the Jaynes-Cummings 
lattice model. Therefore we introduced wave vector and 
filling dependent polariton particle creation and hole cre- 
ation operators, which are linear combinations of photon 
and two- level excitation creation operators. We evalu- 
ated spectral functions and densities of states based on 
the polariton quasiparticles and analyzed the weights of 
their constituents. We have seen that the polariton oper- 
ators are nontrivial combinations of photon and two-level 
system operators, which depend on the wave vector, the 
quasiparticle band, and the filling, or rather the Mott 
lobe. On top of that, the polariton operators of parti- 
cle and hole type are not adjoint operators. It is there- 
fore not possible to describe the JCL model by a simple 
single-band polariton model. 



Appendix A: Properties of the VCA Green's 
function 



Here, we will derive some properties of the VCA 
Green's function for bosons. To begin with, we define 
new operators d\ ~ with ~ = ~ and dl ~ = a~^~ , 

J,r,k l,^,k r,k 2,r,k r,k 

where r stands for the site number within the clusters 
and k is the wave vector of the first Brillouin zone of 
the superlattice. The VCA Green's function will still 
be diagonal in the latter index due to the periodicity of 
the clusters. The spectral representation of the cluster 
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Green's function 

can be written in the compact form using the so-cahed 
Q-matrices2^i2i 

G' = QD:,SQt . 

Here, Q is a M x matrix, where M is twice the number 
of cluster sites (the factor 2 stems from the two species 
of operators) and K = Kp + K^^ where Kp and Kh is 
the dimension of the Hilbert space for Np + 1 and Np — 1 
particles, respectively. The Q-matrix is defined as follows 




(V^ol <^/,r,k for V <K^ 



The diagonal matrix D'^ = diag(cj — contains the 

individual poles uo'^ of the cluster and S = diag(sjy) is a 
diagonal sign matrix with s^y = +1 for particle excitations 
iuo'^ > 0) and Sj^ = —1 for hole excitations {uj^ < 0). 
The VGA Green's function in Q-matrix representation 
for bosons^ reads 

G(k,cj) = QXDc.X-iSQ^ , (Al) 

where D^^ = diag(a; — cjj,)"^ is the diagonal matrix of the 
individual poles at the VGA energies. These energies and 
the corresponding eigenvector matrix X are determined 
via the generalized eigenvalue problem 

(diag(8,c^:,) + QtVQ)X = SXA, 

^ V ' 

= M 

where V = Ho — HqIs the difference of the matrices of 
the single-particle part of the Hamiltonian for the original 
and the reference system (i.e., the cluster). 

A general feature of such eigenvalue equations for Her- 
mitian matrices M is that both X"*" M X and X"*" S X = 
diag(Ai:~^) are diagonalized, but X is not unitary. We 
can exploit this fact as follows 

G(k, a;) = QXD^ X"^ S (X"^)^ X^ . 



The matrix = (X''" SX) ^ is diagonal as we just dis- 
cussed and can be combined with D^; resulting in 

G(k, cj) = QXD, (QX)t 

CO — CO I, 

Moreover, the pole strengths tvi/ are real since 

K;i = (xtsx),, 



When the VGA parameters are determined consistently, 
the stability of the A/p-particle system requires that the 
sign of 1^1, coincides with the sign of the excitation ener- 
gies cjjy, like in the exact spectral representation. 

So far the Green's function still depends on the in- 
tra cluster indices r, s. The purely k-dependent Green's 
function is commonly obtained by Green's function- 
periodization i^^i^^ Invoking the periodization prescrip- 
tion yields the Green's function matrix merely in the 
indices /, J for the two particle species 

G(k, uj) = QXDc^XtQt (A2) 

with Q,^, = l^e-^^^'^Q/,,;,, (A3) 

r 

see also Eq. (f23|) . Equation (|A2[) corresponds to the spec- 
tral representation of the exact Green's function and it 
allows to extract the VGA approximation of the many- 
body eigenstates of the infinite system, which are obvi- 
ously a linear combination of the cluster eigenstates for 
both, particle and hole excitations. 

As described in the text we need the integrated spectral 
density, i. e., 

A/j(k, Qa) = [ (--Im G/j(k, u + if])) duj 

We readily recognize, that the integrated spectral den- 
sity is either positive or negative definite, depending on 
whether the quasiparticle under consideration is of par- 
ticle or hole type. Equivalent ly, in the original represen- 
tation 

A,j(k, n^) = ^ (Q x),,,(x-is Qt),, J . 

For the polariton discussion it is convenient to suppress 
the minus sign arising in the hole case and we define the 
strictly positive integrated spectral densities as 

Aij{k, = |^,j(k, . (A4) 

Appendix B: Solution of the single-site problem 

For zero-hopping strength t = the JGL model can be 
solved exactly, as it reduces to a single-site problem, i. e., 
to the JG model. Including the chemical potential yields 
the single-site Hamiltonian 

^JCL ^ ^JC _ ^ (^t ^ + ^+^-) ^ (Bl) 

where we dropped the site index i. It can be evaluated 
with respect to the bare states |np, 5), where Up is the 
number of photons and s G {|, t}- Next, we sketch 
the most important steps for solving the single-site JGL 
model. A detailed discussion can be found for example 
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in Refs. [25| or [39|. As the JC Hamiltonian conserves the Plugging Eq. (j20aj) in Eq. ()Cip yields 



particle number the Hamiltonian H^^^ is block diago- 
nal. Each block corresponds to a certain particle number 
n and thus we use the bare states \n — 1, t) and |n, 1) to 
evaluate the block, which yields 



{n — l)ujc -\- e — fin ^Jn 

y/n nujc — jJ^n 



(B2) 



when using as denoted in Sec. [Ill the coupling g as unit 
of energy. The eigenvalues of the block are 

E\n,a) = nujc - y -^aq{n) - /an , (B3) 



where a G { — , +} and q{n) = y (A/2) + n. For a 
certain particle number n the energy ^|n,-) is always 
smaller than £^|n,+) and thus ^|n,-) is the ground state 
energy of the sector with n particles, i.e., of the block 
Bn. The eigenvectors |n, a) of the matrix Bn are termed 
dressed states and are given by 



|n, a) = Una \n - 1, t) + ^na i) , 



(B4) 



where n > 0, '^n+) = (sin^(n), cos^(n)) and 

{un-, Vn-) = (cos 6>(n), — sin^(n)) with the following re- 
lations sin 6>(n) = \/ {q{n) — A/2)/2q{n) and cos 6{n) = 
V (^(^) + ^/2)/2g(n). An exception is the bare state 
|0, ^), which forms a 1 x 1 block of zero particles and has 
the eigenvalue £^|o,^) = 0, independently of the detuning 
A. According to the notation used in Eq. (|B4[) . we denote 
this state as |0, — ). In order to obtain the phase bound- 
ary for zero hopping between two adjacent Mott lobes, 
the energies ^|n,-) and E\n-\-i^-) have to be set equal. 
The energies of the states |m, — ) are used, as the phase 
diagram is evaluated for the ground state. The compar- 
ison of the energies yields (/i — ujc) = q{n) — q{n + 1) for 
the location of the phase boundary at zero hopping. 



Appendix C: First-order degenerate perturbation 
theory 

In this appendix we evaluate the results of first-order 
degenerate perturbation theory for the single-particle 
and single-hole excitation bands of the JCL model. To 
apply first-order degenerate perturbation theory the ma- 
trix elements of the perturbation Hi = ^ tij a] aj , 
where tij is the hopping matrix, have to be evaluated 
with respect to the degenerate states and |^^'^), 

see Eq. (|2Q]) . As the hopping term Hi does not change 
the total particle number and does not effect the excita- 
tion a, the following two matrices have to be evaluated; 
one for single-particle excitations 



and one for single-hole excitations 



(Cl) 



N 



(MP«' =(g)(n, -l{n + l,a\iJ2tv4aj 



N 



\n+l,a),Q^\n,-)^, . (C3) 

u' = l 

Due to the orthogonality of the eigenvectors of sectors 
with different particle number, the conditions i = / and 
j = V hold, which reduce the matrix elements to 

{^p)w = Ui' {n, -|^, (n + 1, a\i aj ai' \n + 1, a)^, |n, -)^ 
= tw\{n^l,a\a^\n,-)\K (C4) 

In the second step, we dropped the site index as the ex- 
pectation value does not depend on the specific lattice 
site. The corrected matrix elements are thus the old ones 
with renormalized hopping strength 

-i'^ = -t\{n+l,a\a^\n, -) 1^ 



Analogously, one obtains 

(M^)„. =t,./|(n-l, a\a\n, -) 1^ 



(C6) 



for the matrix elements defined in Eq. (|C2|) . From that 
the renormalized hopping strength for single-hole excita- 
tions is evaluated as 

-t^ = -t\^/n - lUn-lQ, Un- VnVn-lQ.Vn-\^ ' (C7) 

The eigenvalues of the matrices M^^^ are the first-order 
corrections and thus the corrected energies f|n±i, a)(^) of 
the one-dimensional JCL model are given by 

^|n+i, oc) {k) = £^|n+i, a) -^ip COS k and (C8a) 
£\n-i, a) (k) = E\n-i, a) -^ih COS k , (C8b) 

respectively, where /c is a wave vector of the first Brillouin 
zone. Within first-order degenerate perturbation theory 
we obtain 

= ^|n+l, a) {k) - E\n, -) 

= {uc - fi) ^ a q{n + 1) + q{n) -2ip cos k (C9) 
for the single-particle excitation band and 

^h,l = ^|n, -) - £\n-l, a) (^) 

= {uJc- 1^) - a q{n - 1) - q{n) + 2 ?^ cos A: . (CIO) 



(M-V^(^^'^|^i|^^'^). 



(C2) for the singe-hole excitation band. 
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Appendix D: Polariton operators in the single-site 
limit 



With the definition Xq, = (^Qi^c^, ^Q2,a)^ integrated 
spectral density for the particle part can be expressed as 



In this appendix, we want to analyze the polaritonic 
feature in the single-site limit for zero detuning. 

In the single-site limit it is exactly possible to construct 
a polariton operator which, applied to the many-body 
eigenstate |n, — ), generates the eigenstates |n±l,a). 
The polaritonic weights follow from 



(/?at+7a+)|n,-) 



PI 



V2 

I 

= |n + l,a) 
1 



\nA)-p\l'^\n^lA) 



(Dl) 



Here, we explicitly include the filling n as index. So the 
relative weights are (V^ + 1 + V^)~"^ upper band 

(a = +) and — (^/rH^+y^) lower band (a = — ). 

This means that in the lower band the weights have oppo- 
site sign and the polaritons are of predominant photonic 
character, while the opposite applies to the upper band. 
The modulus of relative weight is just the inverse, i.e., 

l/3+Mi = i7-/ri. 

Next, we study the hole case for n > 1 
(/3a + 7cr") |n, -) = 



P\l^r- \n - 2,t) + ^ 7 1^ - 1.^) 



V2 



= |n — 1, a) 

^In _ 1 



\/n + ay/n - 1 ' 



(D2) 



which is positive for both bands a = ±. Again we have 
the reciprocal property /3+/7+ = 7~//3~ and the lower 
band has predominantly photonic character, while the 
opposite is the case in the upper band. 

Now we want to scrutinize the generalized eigenvalue 
problem of the Green's function. The single-site Green's 
function reads 



E 



'Qi,^ = {n^ha\d\\n,-y 



1, a\ dj In, - 



where we introduced the operators di = a and d2 = cr 
For the single-particle term we obtain 



The overlap matrix Sp = {djd^j) is readily obtained by 
the spectral theorem 

Sp = x+ x+ + x_ x_ , 

and the generalized eigenvalue problem for the polariton 
weights according to Eq. ([15]) reads 



(1 A) Xq, Xq, Zq, — Ax_Q, x_Q, Zq, . 

The eigenvalues are zero and one. For the polariton 
weights we are interested in the latter. The correspond- 
ing eigenvector is simply given by the vector orthogonal 
to X_f> 



2 \y/n + ay/n + 1 / 
With that one obtains for the ratio of the weights 

0"^ 1 

7p,n + a^/nTT ' 

which is in agreement with the exact result in Eq. ()Dip . 
Now we address the hole case, again for n > 1, 

1 



Qi^a = {n-l,a\a |n, -) -(V^ - 1 - a^/rl) 

^Q2,a = (n + 1, a\ a~ |n, -) = . 

We proceed as in the particle case with the definition 
of Xq, ^ = (^Qi^Q, ^Q2,a)' The remaining steps are the 
same as before and we end up with 



PI 



a 



7^^^ \/n — 1 + a^/n yjn + OL\/n — 1 

which is also in agreement with the exact result, see 
Eq. ()D2p . So we see that the determination of the po- 
laritonic weight via the generalized eigenvalue problem is 
reasonable. In the single-site limit, the exact many-body 
eigenstates |n =b 1, a) can be generated correctly by suit- 
able polariton operators acting on the state |n, — ). The 
operators are, however, not universal, they depend on n 
and in the lattice case even on k. On top of that, the 
polariton creation operator for holes is not the adjoint 
of the corresponding polariton creation operator of the 
particle type, or in other words its annihilation operator. 
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